{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Imports\n",
    "\n",
    "# utility modules\n",
    "import glob\n",
    "import os\n",
    "import sys\n",
    "import re\n",
    "\n",
    "# the usual suspects:\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# specialty modules\n",
    "import h5py\n",
    "import pyproj\n",
    "\n",
    "# run matplotlib in 'widget' mode\n",
    "%matplotlib widget\n",
    "%load_ext autoreload\n",
    "%autoreload 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "sys.path.append(os.path.join(os.getcwd(), '../..'))\n",
    "from readers.read_HDF5_ATL03 import read_HDF5_ATL03\n",
    "from readers.get_ATL03_x_atc import get_ATL03_x_atc"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def atl06_to_dict(filename, beam, field_dict=None, index=None, epsg=None):\n",
    "    \"\"\"\n",
    "        Read selected datasets from an ATL06 file\n",
    "\n",
    "        Input arguments:\n",
    "            filename: ATl06 file to read\n",
    "            beam: a string specifying which beam is to be read (ex: gt1l, gt1r, gt2l, etc)\n",
    "            field_dict: A dictinary describing the fields to be read\n",
    "                    keys give the group names to be read, \n",
    "                    entries are lists of datasets within the groups\n",
    "            index: which entries in each field to read\n",
    "            epsg: an EPSG code specifying a projection (see www.epsg.org).  Good choices are:\n",
    "                for Greenland, 3413 (polar stereographic projection, with Greenland along the Y axis)\n",
    "                for Antarctica, 3031 (polar stereographic projection, centered on the Pouth Pole)\n",
    "        Output argument:\n",
    "            D6: dictionary containing ATL06 data.  Each dataset in \n",
    "                dataset_dict has its own entry in D6.  Each dataset \n",
    "                in D6 contains a numpy array containing the \n",
    "                data\n",
    "    \"\"\"\n",
    "    if field_dict is None:\n",
    "        field_dict={None:['latitude','longitude','h_li', 'atl06_quality_summary'],\\\n",
    "                    'ground_track':['x_atc','y_atc'],\\\n",
    "                    'fit_statistics':['dh_fit_dx', 'dh_fit_dy']}\n",
    "    D={}\n",
    "    file_re=re.compile('ATL06_(?P<date>\\d+)_(?P<rgt>\\d\\d\\d\\d)(?P<cycle>\\d\\d)(?P<region>\\d\\d)_(?P<release>\\d\\d\\d)_(?P<version>\\d\\d).h5')\n",
    "    with h5py.File(filename,'r') as h5f:\n",
    "        for key in field_dict:\n",
    "            for ds in field_dict[key]:\n",
    "                if key is not None:\n",
    "                    ds_name=beam+'/land_ice_segments/'+key+'/'+ds\n",
    "                else:\n",
    "                    ds_name=beam+'/land_ice_segments/'+ds\n",
    "                if index is not None:\n",
    "                    D[ds]=np.array(h5f[ds_name][index])\n",
    "                else:\n",
    "                    D[ds]=np.array(h5f[ds_name])\n",
    "                if '_FillValue' in h5f[ds_name].attrs:\n",
    "                    bad_vals=D[ds]==h5f[ds_name].attrs['_FillValue']\n",
    "                    D[ds]=D[ds].astype(float)\n",
    "                    D[ds][bad_vals]=np.NaN\n",
    "    if epsg is not None:\n",
    "        xy=np.array(pyproj.proj.Proj(epsg)(D['longitude'], D['latitude']))\n",
    "        D['x']=xy[0,:].reshape(D['latitude'].shape)\n",
    "        D['y']=xy[1,:].reshape(D['latitude'].shape)\n",
    "    temp=file_re.search(filename)\n",
    "    D['rgt']=int(temp['rgt'])\n",
    "    D['cycle']=int(temp['cycle'])\n",
    "    D['beam']=beam\n",
    "    return D"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "D_2l={}\n",
    "D_2r={}\n",
    "\n",
    "data_root='/srv/shared/surface_velocity/'\n",
    "field_dict={None:['delta_time','latitude','longitude','h_li', 'atl06_quality_summary'],\\\n",
    "                    'ground_track':['x_atc','y_atc'],\\\n",
    "                    'fit_statistics':['dh_fit_dx', 'dh_fit_dy']}\n",
    "\n",
    "# specify the rgt here:\n",
    "rgt=\"1092\"\n",
    "# iterate over the repeat cycles\n",
    "for cycle in ['03','04','05','06','07']:\n",
    "    for filename in glob.glob(os.path.join(data_root, 'FIS_ATL06_small', f'*ATL06_*_{rgt}{cycle}*_003*.h5')):\n",
    "        try:\n",
    "            # read the left-beam data\n",
    "            D_2l[filename]=atl06_to_dict(filename,'/gt2l', field_dict=field_dict, index=None, epsg=3031)\n",
    "            # read the right-beam data\n",
    "            D_2r[filename]=atl06_to_dict(filename,'/gt2r', field_dict=field_dict, index=None, epsg=3031)\n",
    "            # plot the locations in the previous plot\n",
    "            #map_ax.plot(D_2r[filename]['x'], D_2r[filename]['y'],'k');  \n",
    "            #map_ax.plot(D_2l[filename]['x'], D_2l[filename]['y'],'k');\n",
    "        except Exception as e:\n",
    "            print(f'filename={filename}, exception={e}')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "b6cbb7b6cdd74e52a18393fd0491ab8e",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7fea8b0d1290>]"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "for i, key in enumerate(D_2l.keys()):\n",
    "    if i == 0:\n",
    "        D1 = D_2l[key]\n",
    "            \n",
    "for i, key in enumerate(D_2l.keys()):\n",
    "    if i == 2:\n",
    "        D2 = D_2l[key]\n",
    "\n",
    "#for filename, Di in D_2l.items():\n",
    "#Plot only points that have ATL06_quality_summary==0 (good points)\n",
    "ind1=D1['atl06_quality_summary']==0\n",
    "ind2=D2['atl06_quality_summary']==0\n",
    "#define the heights of the segment endpoints.  Leave a row of NaNs so that the endpoints don't get joined\n",
    "h1 = D1['h_li'][ind1]\n",
    "x1 = D1['x_atc'][ind1]\n",
    "h2 = D2['h_li'][ind2]\n",
    "x2 = D2['x_atc'][ind2]\n",
    "\n",
    "plt.figure()\n",
    "plt.plot(x1,h1)\n",
    "plt.plot(x2,h2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "a233258343e94ffcab153f61ef0d5ae1",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "idx = ~np.isnan(h1) & ~np.isnan(x1) & ~np.isnan(h2) & ~np.isnan(x2)\n",
    "H1 = h1[idx]\n",
    "X1 = x1[idx]\n",
    "Slope1 = np.gradient(H1,X1)\n",
    "H2 = h2[idx]\n",
    "X2 = x2[idx]\n",
    "Slope2 = np.gradient(H2,X2)\n",
    "\n",
    "start, end = 500,1000\n",
    "\n",
    "plt.figure();\n",
    "plt.subplot(211)\n",
    "plt.plot(X1[start:end],H1[start:end], label=\"cycle=3\")\n",
    "plt.plot(X2[start:end],H2[start:end], label=\"cycle=5\")\n",
    "plt.xlabel('x_atc')\n",
    "plt.ylabel('elevation');\n",
    "plt.legend()\n",
    "\n",
    "plt.subplot(212)\n",
    "plt.plot(X1[start:end],Slope1[start:end],ms=1, label=\"cycle=3\")\n",
    "plt.plot(X2[start:end],Slope2[start:end],ms=1, label=\"cycle=5\")\n",
    "plt.xlabel('x_atc')\n",
    "plt.ylabel('slope');\n",
    "plt.legend()\n",
    "\n",
    "\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "18010e94870c4faa86d6d7d7661c7268",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "velocity of the best fit:  1.095523313877746 m/d\n",
      "velocity of the best fit:  399.8660095653773 m/a\n"
     ]
    }
   ],
   "source": [
    "start = 500\n",
    "L = 500\n",
    "mdpt = int(start+L/2)\n",
    "xctr = X1[mdpt]\n",
    "sec_offset = abs(D1['delta_time'][mdpt]-D2['delta_time'][mdpt])\n",
    "day_offset = sec_offset/(3600.*24.)\n",
    "x_step = np.nanmean(np.gradient(X1[start:start+L]))\n",
    "\n",
    "idx_offsets = np.arange(-100,100)\n",
    "cs = np.empty_like(idx_offsets).astype(float)\n",
    "for j,offset in enumerate(idx_offsets):\n",
    "    cs[j] = np.corrcoef(Slope1[start:start+L],Slope2[start+offset:start+L+offset])[1,0]\n",
    "\n",
    "idx_best = idx_offsets[np.argmax(cs)]\n",
    "N_offset = abs(idx_best)\n",
    "dx_offset = x_step*N_offset\n",
    "vel = dx_offset/day_offset\n",
    "\n",
    "plt.figure();\n",
    "plt.subplot(211)\n",
    "plt.plot(X1[start:start+L],H1[start:start+L], label=\"cycle=3\")\n",
    "plt.plot(X2[start:start+L],H2[start+idx_best:start+L+idx_best], label=\"cycle=5\")\n",
    "plt.xlabel('x_atc')\n",
    "plt.ylabel('elevation');\n",
    "plt.legend()\n",
    "\n",
    "plt.subplot(212)\n",
    "plt.plot(X1[start:start+L],Slope1[start:start+L],ms=1, label=\"cycle=3\")\n",
    "plt.plot(X2[start:start+L],Slope2[start+idx_best:start+L+idx_best],ms=1, label=\"cycle=5\")\n",
    "plt.xlabel('x_atc')\n",
    "plt.ylabel('slope');\n",
    "plt.legend()\n",
    "\n",
    "plt.tight_layout()\n",
    "\n",
    "print('velocity of the best fit: ' , dx_offset/day_offset, 'm/d')\n",
    "print('velocity of the best fit: ' , (dx_offset/day_offset)*365, 'm/a')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "starts = np.arange(500,len(X1)-600,100)\n",
    "L = 500\n",
    "\n",
    "idx_offsets = np.arange(-100,100)\n",
    "vels = np.empty_like(starts).astype(float)\n",
    "xctrs = np.empty_like(starts).astype(float)\n",
    "\n",
    "for i,start in enumerate(starts):\n",
    "    mdpt = int((start+end)/2)\n",
    "    xctrs[i] = X1[mdpt]\n",
    "    sec_offset = abs(D1['delta_time'][mdpt]-D2['delta_time'][mdpt])\n",
    "    day_offset = sec_offset/(3600.*24.)\n",
    "    x_step = np.nanmean(np.gradient(X1[start:start+L]))\n",
    "\n",
    "    cs = np.empty_like(idx_offsets).astype(float)\n",
    "    for j,offset in enumerate(idx_offsets):\n",
    "        cs[j] = np.corrcoef(Slope1[start:start+L],Slope2[start+offset:start+L+offset])[1,0]\n",
    "    \n",
    "    idx_best = idx_offsets[np.argmax(cs)]\n",
    "    N_offset = abs(idx_best)\n",
    "    dx_offset = x_step*N_offset\n",
    "    vels[i] = dx_offset/day_offset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "40b8bb2cc0794e1da2b91a35fce2199c",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure()\n",
    "plt.subplot(211)\n",
    "plt.plot(xctrs,vels*365.)\n",
    "plt.ylabel('Velocity (m/yr)')\n",
    "plt.subplot(212)\n",
    "plt.plot(x1,h1)\n",
    "plt.plot(x2,h2)\n",
    "plt.ylabel('Elevation (m)')\n",
    "\n",
    "plt.tight_layout()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
